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Abstract 

Hard-sphere fluids confined between paraUel plates a distance D apart are studied for a wide range of 
packing fractions, including also the onset of crystallization, applying Monte Carlo simulation techniques 
and density functional theory. The walls repel the hard spheres (of diameter a) with a Weeks-Chandler- 
Andersen (WCA) potential Vwca{z) = 4e[((7^/z)^^ — (cr^/z)^ + 1/4], with range aw = o l2. We vary the 
strength e over a wide range and the case of simple hard walls is also treated for comparison. By the 
variation of e one can change both the surface excess packing fraction and the wall-fluid (7^;/) and wall- 
crystal (7«,c) surface free energies. Several different methods to extract 7^/ and from Monte Carlo 
(MC) simulations are implemented, and their accuracy and efficiency is comparatively discussed. The 
density functional theory (DFT) using Fundamental Measure functionals is found to be quantitatively 
accurate over a wide range of packing fractions; small deviations between DFT and MC near the fluid 
to crystal transition need to be studied further. Our results on density profiles near soft walls could be 
useful to interpret corresponding experiments with suitable colloidal dispersions. 
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I. INTRODUCTION 



Recently it has been demonstrated that for colloidal suspensions t 
tunable from hard spheres to soft repulsion jl-3| or weak attraction 



the structure of fluid-crystal 



0|io| 



le effective interactions are 



6[, and at the same time 



and fluid-wall interfaces can be analyzed in arbitrary detail, 



10|. Correspondingly, there is a 



e.g. by visualizing the packing of particles in these interfaces 
great interest in model studies pertinent to such systems. However, most work has focused on 
the archetypical hard sphere fluid 
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14| . confined by hard walls 



15 



28 



heterogeneous crystal nucleation at hard walls 



23 



30|. With respect to 



24| . this system is difficult to understand, since 



there is evidence that complete wetting of the wall by the crystal occurs, when the fluid packing 
fraction approaches the fluid-crystal phase boundary in the bulk [24 1. 



Now it is well known that the interaction between colloidal particles and walls can also be 
manipulated, by suitable coatings of the latter, e.g. via a grafted polymeric layer (using the 
grafting density and chain length of these polymers, under good solvent conditions, as control 



parameters 



31 



32|). Thus, in the present work we explore a model where colloidal particles that 



have an effective hard-sphere interaction in the bulk experience a soft repulsion from confining 



33| 



walls, describing this repulsion for the sake of simplicity by the Weeks-Chandler-Andersen 
potential. We show that such a short-range repulsion has only small effects on the structure of 
the fluid near the wall, but nevertheless affects the wall-fluid interface tension 'jyjf significantly. 
Both Monte Carlo methods and density functional calculations are used. 



In Sec. 2, the model is introduced, and several Monte Carlo methods to extract 7,„/- a.re briefly 



described. Since the judgment of accuracy for such methods is somewhat subtle (28[, we are 
interested in comparing estimates from several rather different approaches, to avoid misleading 
conclusions. In Sec. 3, we present our results on density profiles, while Sec. 4 describes our 
results for the dependence of 7^,/ on packing fraction. First preliminary results on the interfacial 
tension between the crystalline phase and the confining wall are presented in Sec. 5, while Sec. 6 
summarizes our results and discusses possible applications to experiments. The density functional 
methods are briefly explained in an Appendix. 



II. MODEL AND SUMMARY OF MONTE CARLO METHODS FOR THE ESTIMA- 
TION OF WALL FREE ENERGIES 



The simulated model is the simple fluid of hard particles of diameter a, in the geometry of an 
Lx Lx D system, confined between two parallel walls located at z = and a.t z = D. In the x and 
y directions, periodic boundary conditions are applied throughout. The particle- wall interaction 
contains either a hard wall type interaction 

Vhw{z) = oo for z < cr/2 and for z > D — a/2 (1) 



or a soft repulsion of the Weeks- Chandler- Andersen [33| type 



Vwca{z) = 4e [{a^/zf^ - {a.^/zf + i] for < z < a^2'l^ 

= 4e [{a^/iD - z)^ - i^n^/iD - z)f + \] for {D - a^2^l^) <z<D (2) 
= otherwise 
In Eq. [2] we choose = u /2, while the parameter e that controls the strength of this additional 
soft repulsion is varied in the range < e < 4 (choosing units such that ks = I and absolute 
temperature T = 1). Note that for e = Eq. [2] also becomes a hard-core potential Vhc{z) = oo 
for 2 < and z > D, respectively. 

In the literature the confinement of hard spheres between hard walls, i.e. the case where only 



Vhw{z) is present, has already been extensively studied |l5l-l28l. l30| . while we are not aware of 
any work using Vwca{z) instead. The advantage of the choice Eq. [2] from the theoretical point of 
view, is that e is a convenient control parameter: varying e the wall-fluid interfacial tension 7^^- 
as well as the wall-crystal interfacial tension can be modified. Note that the direct effect of 
Vwca{z) is zero in the range au,2^/^ < z < D — a^2^/^: thus, when D is very large, we expect 
that the structure of the hard sphere fluid in the center of the slit (very far from both walls) is 
identical to a corresponding hard sphere fluid in the absence of confining walls (applying periodic 
boundary conditions also in the ^-direction). 

We stress that the WCA form of the potential in Eq. [2] is only chosen for the sake of computa- 
tional convenience. Having the application to colloidal dispersions in mind, one might expect that 
the colloidal particles carry a weak electrical charge, but the Coulomb interactions are strongly 
screened by counterions in the solution. Assuming also some effective charges at the walls, a po- 
tential like C exp(— Kz) might seem a physically more natural choice (with a screening length 



of the order of a/ 10 [3| or even smaller; note that the constant C could be positive or negative). 
However, when flexible polymers are adsorbed (or end-grafted) at the walls, the chain length 
and grafting density cr^ provide additional parameters of a repulsive potential due to the dangling 
chain ends out in the solution, if there is no adsorption of the polymers on the colloidal particles. 
Thus, the actual potential between colloidal particles and confining walls is clearly non-universal, 
it depends on system preparation and can be fairly complicated due to a superposition of several 
mechanisms. Since we do not attempt to model any specific system, we take Eqs. [Tpi as a generic 
model. 

For simulations in the standard canonical (constant volume) ensemble, the standard Monte 



Carlo algorithm 3J] with local single particle moves is implemented, choosing particles at random 
and attempting to move their center of mass to a new position. Of course, moves are accepted 
only if they respect the excluded volume between the particles. For the system with walls, the 
Metropolis criterion needs to be tested if either the old or the new position of the particle is within 
the range of the wall potential, Eq. [2l At this point, the advantage of choosing a potential that is 
strictly zero for a broad range of z (as specified above) clearly becomes apparent. 

The observables of interest (for simulations in the canonical ensemble) are the normal pressure 
Pjv and the local tangential pressure Pt{z) and the corresponding number density profile p{z), 

D 

choosing the average particle density p = f p{z)dz/D = N/{L'^D) or the corresponding packing 



fraction 

ri = (vraV6)p (3) 
as the input parameters that we vary in our simulation. 

Note that due to wall effects on the hard sphere fluid we expect an approach to the bulk density 



Pb(P/v,T) as D — J- oo as follows [35| 



p = p,{Pn, T) + 2pjD, D^oo (4) 

where the surface excess density ps (and associated surface excess packing fraction 77^) are 
formally defined for a semi-infinite system as 

00 

= y [P(^) - Pb]dz , r]s = Ps'rr/6 . (5) 



In a film of finite thickness -D, an analog of Eq. [S]can be used if p{z) has settled down to pb 
already for values of z that are clearly smaller than D/2: then the upper limit oo in Eq. [5] can be 
replaced by D/2 with negligible error. In this limit, the two walls can be considered as strictly 



35| simply related to the 



non-interacting, and then the wall-fluid interfacial tension is also 

D 

difference between Pjy and the average tangential pressure, Pt = J dzPj'{z) / 



7«./ = {Pn - Pt)D/2. (6) 

Note, however, that the situation is more subtle for a crystal confined between two walls, since 
the long range crystalline order in the crystal is not necessarily commensurate with the chosen 
distance D and hence the long-range elastic distortion of the crystal that will in general result 
invalidates the above statement that the effects of the two walls add independently. But, for fluid 
systems Eq. [H] is useful if D is large enough. 



As is we 
expression 



1 known, the standard "mechanical" approach to calculate the pressure from the virial 



35 



371 ] cannot be straightforwardly applied for systems with hard-core interactions. In 
order to apply Eq. [6l we thus follow the approach of de Miguel and Jackson We here recall 
only briefly the most salient features. For a bulk hard sphere fluid the number of pairs with a 
relative distance in the range from a to cr + Ar is sampled, n(Ar), and one estimates the derivative 



of this function for Ar — )• 0, a = a'd{n{Ar))/dr, and uses the formula 27|] 



P/ipkBT) = l + a/i3N) (7) 

to obtain the (average) pressure of a bulk hard sphere fluid at given density p = N/V. Alterna- 
tively, one can consider virtual volume changes by a factor ^ and compute the probability PnoviO 
that there are no molecular pair overlaps when the volume is decreased from V to V = V{1 — 
For small ^ one can show that PnoviO — 6xp(— 6^), where 6 > is related to the pressure by a 
relation similar to Eq. [7] 

P/{pkBT) = 1 + b/N, (8) 

and one can numerically verify that both routes based on Eqs. [71 [8] work in practice, and 
agree within their statistical errors. The method of Eq. [8] now can be straightforwardly extended 
to sample Pn and Pt separately: one considers volume changes that are due to reducing the 



distance from D to D' = D{1 — ^) keeping the lateral distance L constant to obtain P^, while 
L' = L(l — at fixed D is used to obtain Pt {27] . 

When we vary the strength e of the WCA {Eq. [2I[ wall potential for fixed total particle number 
and fixed linear dimensions L and D, the change of ps caused by the variation of e necessarily 
cause a change of pb (and hence P/v), since in the canonical ensemble the total density is strictly 
constant. However, this effect is clearly undesirable: we want to vary e and ps but keep the bulk 
conditions unchanged! Hence it would be preferable to vary e and keep P/v constant, rather than 
keeping the volume constant. But we do wish to keep D constant as well. At first sight, one 
might conclude that these constraints are impossible to realize, since P/v and D are a pair of 



thermodynamically conjugate variables. However, Varnik [38| has devised an iterative method, 
where only the area A = L"^ rather than the whole volume V = L^D is allowed to fiuctuate, as it 
would happen in an NPT ensemble. Applying this method (for details, see jssi) one can realize a 
NP]S[DT ensemble, and this ensemble is indeed useful to implement the variation of e. However, 
due to the larger computational effort of this method our simulations were done in the canonical 
NLDT ensemble. 

In some cases of interest it suffices to compute differences A7 = 7io/(e) — 7«,/(eo) only, rather 
than the values 7u,/(e), 7«,/(eo) individually. E.g., for cq = our model reduces to the case 



of hard walls only, Eq. [H which has been studied extensive 



rather precise values of 7^/(00) are already available 2m |25|. Such differences A7 can, at least 



y in the literature 



15 



-28 



30|, and 



in principle, be found from a thermodynamic integration method based on linear response theory. 
We note that the thermodynamic potential can be written as (for A^ — )■ 00, D — )■ 00) 



G(P^, A^, P), T) = -ksT In j dX exp | - ^biX) - (3PnL^D - 

D 

-Pel' j p{z,X)V{^cA{^)dz^ (9) 


Here prefactors of the partition function that are unimportant for the following argument are 
already omitted. /3 = l/ksT and X stands for a point in the configuration space of the system 
(i.e., X is just the set of coordinates of all the center of masses of the hard spheres). By 'Hb{X) we 
denote the interaction among the hard spheres (i.e., exp[— /^^^^(X)] = if any pair of hard spheres 
overlaps). The interaction with the WCA-potential has been written out explicitly, denoting 



Vwca{z) = ^V{nf(jj^{z), and defining p{z,X) as the particle density in the infinitesimal interval 
[z, z + dz\ for the configuration X. 



We now consider the derivative of G with respect to e, to find (see also 39|] for a related 



treatment of a binary Lennard- Jones mixture) that this derivative just can be interpreted as the 

dt 

D/2 



sum of L^(^^$^) for the two walls (which are identical). Hence 



97. 



wf 



de 



PnDT 





{p{z,X)),V^cA(^)dz, (10) 



where the notation p{z, e) = {p{z,X))e is used to emphasize that the statistical average {■ ■ ■)e 
is carried out in an ensemble where a wall potential Vwca{z) = ^V^^YCAi'^) nonzero only for 
< 2; < 2^/'^a^ = 0.5 ■ 2^/^ ^ 0.561, for our choice a^^ = a/2 = 1/2. From Eq. [lO] we realize 
that any change of '-/wf due to the variation of e can only be due to the fact that the product 
p{z,e)V^Q^{z) changes when e is varied. Now differences A7 can be computed from 



D/2 

df=! df' 

/ PnDT 

eo 



^7 = / ( %1 de= [ de' I p{z,e')V{,aAi^)dz (11) 



For use in actual computations the application of Eq. [TT] is subtle since one needs to record 
p{z, e') in the range < 2 < 2^/^(T^ with very high precision. At this point, we draw attention to 



another version of thermodynamic integration (termed "Gibbs-Cahn integration" 30|) which can 
also be implemented if only hard walls are present (Eq. [1]): there one uses ps {Eq. [5l[ to study the 
variation of 7^/ with the bulk density pb{PN,T) of the system (see Eq. [H] below). 

However, in the present work we rather use another variant of thermodynamic integration, 
which is briefly characterized below. This method (which we refer to as "ensemble switch method" ) 
is a variant of the method used by Heni and Lowen [l^, where one gradually switches between 
a system without walls, applying periodic boundary conditions throughout, described by Hamil- 
tonian 'Hi(X), and a system (with the same particle number and volume) with walls, 'H2{X), 
writing the total Hamiltonian as 

V,{X) = {I - k)V.i{X) + kV.2{X) (12) 

where k G [0, 1] is the parameter that is varied for calculating the free energy difference between 
the systems (1,2). In a simulation k is typically discretized and the system is allowed to move 



from Ki to or with a Metropolis step. The free energy between the two states is given 
by kBT{\nP{i) — In P{i — 1, i + 1)) where P{i) is the relative probability of residin g in state 



i. As this probability varies considerably with k, a variant of Wang-Landau sampling j40|, |4l| is 
employed to eventually simulate each state with equal probability. In this way we can sample free 
energy differences relative to the free energy of the system with periodic boundary conditions as 



a function of k (for technical details see 42|]). For k = 1 the wall free energy then follows from 



7./(p,T)= hm (13) 

with A being the area of the wall. 

Note that also in this method for finite D the density p in the system with walls (k = 1) differs 
from the corresponding system with periodic boundary conditions (k = 0) due to the surface 
excess density. Thus an extrapolation to D — t- oo is necessary. 

III. DENSITY PROFILES OF HARD SPHERE FLUIDS CONFINED BETWEEN WCA 
WALLS 

Figs. [1], [2] show typical data for the density profiles p{z) obtained from our simulations, using a 
box of linear dimensions L = 12.41786, D = 25.61184, and varying the particle number as well 
as the strength e of the WCA potential, Eq. O 

At first sight, the density profiles for the different choices of e look essentially identical; only 
when a magnified picture of the first peak of p{z) adjacent to one of the walls is taken, one sees 
a systematic effect: the larger e, the more remote from the wall the peak occurs, as expected. 
However, for o"„2^/^ < z < D — a^^^^^, i.e. outside the range where the wall potential acts, the 
effect of varying e is negligible. However, for packing fractions t] close to the value rjh cr where 
in the bulk crystallization starts to set in, rjb^cr = 0.492, such as rjb ~ 0.47 or larger, the wall- 
induced oscillations in the density profile ("layering") extend throughout the film (Fig. [2]). This 
observation indicates that the chosen thickness D, as quoted above, is not large enough to allow 
an approach very close to the transition, when one tries to disentangle the effects of the walls (as 
measured by ps or 77^, Eq. [5l respectively) and pb {Eq. H]} or 77^. 

We have compared the values for the normal pressure P/v and corresponding value of rjb as 
function of the nominal packing fraction (77) chosen in our simulations, for a range of values for 



e, the strength of the WCA potential at the walls to literature data 27|, l28|, ISOj. This shows 



that in the chosen range of rjb the hnear dimensions L and D chosen here are large enough to 
allow a meaningful estimation of rjb- Due to the surface excess of the density, there is a systematic 
discrepancy between 77^ (the packing fraction in the center of the thin film) and rj (the total packing 
fraction in the film). 

Fig. [3] shows a plot of rjs (which turns out to be negative for all parameters that were studied) 
versus rjb. Corresponding results from the DFT calculations (see the Appendix for technical details) 
are included. One sees that rjs depends in a nontrivial way on both r/^ and e. It can also be seen 
that for 77buik ~ 0.4 systematic discrepancies between DFT and simulation start to occur, while for 
smaller ?7buik both methods are in excellent agreement. Interestingly, for e = 1 the data are rather 
close to the case where a hardcore potential is used at the walls (data labeled as HW in Fig. [3]). 
The latter case has been studied before by Laird and Davidchack 30|, and the present calculation 
is found to be in excellent agreement with these recent results. This very good agreement is rather 
gratifying, since the latter authors have studied a much larger system {D = 65a, L = 50a) than 
we have used. However, such larger systems are needed very close to the liquid-solid transition, 
due to the extended range of the layering (Fig. [2]). It is also suggestive that the behavior of rjg for 
e — )■ is singular (this limit again corresponds to the hard wall case, but a hard wall at a position 
shifted by a/2). Note that the choice of the square cross section of the box (together with the 
periodic boundary condition) does not lead to noticeable systematic errors. Computations with 
a rectangular Lx x Ly cross section (compatible with a perfect triangular lattice of close-packed 
planes parallel to the walls) have also been made, but the results agree with those that are shown 
within the size of the symbols. 

The distinct effect of the variation of e on 7]^ seen in Fig. [3] can already be taken as an indication 
that a clear effect on the interfacial tension 7^^ can also be expected. To elucidate this point 
further, we present in Figs. HI [5] in more detail the behavior of both p{z) and ep{z,e)V^Qy^{z). 
Recall that the product p{z, e')^iycyi(^) appears in the integral when we relate 7(77?,, eo) and 7(?7b; e) 
by thermodynamic integration {Eq. Illlj . Indeed one can see that the functions p[z,e)VwcA{z) 
does depend on e significantly. 

However, it is also clear from Fig. H] that the use of Eq. [11] for practical computations would be 
difficult, since a very fine resolution of the z-dependence is necessary (while z varies in between 
< z < D = 25.61184, the important intervals contributing to Eq.fTTlhave only a width Az ^ 0.1, 
and the location where these important intervals occur depend on e and are not known precisely 



beforehand). Nevertheless the data of Figs. [31 IH show that varying e does have a pronounced effect 
on both the surface excess density ps (or packing fraction 77^, respectively) and on the function 
Vwca{z)p{z), and hence it is clear that varying e must lead to a change of 7^/- as well. This will 
be explored in the next section. It is also very gratifying that with respect to Vwca{z)p{z), the 
quantity that controls the surface tension 7«,/(e), there is excellent agreement between the MC 
estimates and the DFT calculations for a wide range of packing fractions (0 < rjh < 0.45). Only in 



the immediate vicinity of the freezing transition (0.46 < rji, < rjf, with [13|, |43| rjf = 0.492) slight 
but systematic deviations are apparent in Fig. Hb. For the surface excess density ps, however, 
which is sensitive to the whole profile p{z) and not only to the peaks of p{z) next to the walls, 
deviations between DFT and MC start at smaller r/f, already. We add the caveat, however, that 
close to freezing the finite size effects on the density profile p{z) need to be carefully studied (see 
the discussion of Fig. [2]) but this is left to future work. 

IV. SURFACE FREE ENERGIES OF THE HARD SPHERE MODEL IN THE FLUID 
PHASE 

As a test of our MC procedures, it is useful again to consider the hard wall case {Eq. [TO] } 



first, since this case has been extensively studied in the literature [19[, [27|, [28|, [30|. Fig. [5] gives 
evidence that our methods (based on Eq. [6] or Eq. [121 respectively) are in mutual agreement 
and in agreement with the calculations in the literature, within the statistical errors expected for 
these data. Again the DFT calculation is in excellent agreement over a wide range of packing 
fractions with the simulation results. Only close to the freezing transition small but systematic 
deviations are present, as can be expected from the differences in the surface excess density p^ 
close to freezing (see Fig. [3]). The surface excess density and the surface tension are connected 
through the Gibbs adsorption relation 

dpb oPN{pb) 

where p^ and P/v are the chemical potential and the bulk (normal) pressure, respectively, pertaining 
to the bulk density p^. 

Having asserted that the errors of our calculations are reasonably under control, for the standard 
hard wall case, we turn to the problem of main interest in the present work, namely the variation 
of 7«,/(e) with the strength e of the WCA potential (Figj6]). As we had expected, by changing e we 
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can indeed obtain a variation of 7«,/(e) over a wide range. It is slightly disturbing, however, that 
there seem to be slight but systematic discrepancies between the MC results obtained from Eq. [6] 
and those from the thermodynamic integration method, Eq. [121 this shows that the judgment of 
systematic and statistical errors in these methods is somewhat subtle. However, if we allow for 
statistical errors of the order of three standard deviation rather than one standard deviation, there 
would no longer be any significant discrepancy. Since for most purposes such a moderate accuracy 
in the estimation of 7«,/(e) is good enough, we have not attempted to significantly improve the 
accuracy of our simulations, since this would require a massive investment of computer resources. 
Finally, we note that again the DFT results are very close to the MC data, particularly for 
rjb < 0.45 while closer to the freezing transition small but systematic discrepancies occur again. 
This very good agreement between DFT and simulations for 7^/ is expected from the fact that 
DFT describes the density very accurately close to the walls where Vwca acts. More prominent 
deviations in the density profiles between simulation and DFT are seen near the second peak from 
the wall. DFT does not seem to account for its precise shape near freezing. This deficiency is also 
visible in the "hump" in the second peak of the pair correlation function near freezing which can 



be interpreted as a structural precursor to the freezing transition 



45|. 



V. SOME RESULTS ON THE WALL-CRYSTAL SURFACE TENSION 

As has already been stated earlier, studying the wall-crystal surface free energy is a subtle 
matter, since (i) in general there is always a misfit in a thin film geometry between the distance D 
between the walls, and the lattice spacing a{rib) which depends on the packing fraction in the bulk, 
of course. In addition (ii) the wall-crystal free energy depends on the orientation of the crystal 
axes relative to the walls. In the present context, it is natural to restrict attention to a crystal 
orientation only where the close packed (111) planes at the face-centered cubic crystal structure 
(remember that in the fcc-structure there is an ABCABC... stacking of close-packed planes having 
a perfect triangular crystal structure each) are parallel to the planes forming the walls. Of course, 
it is this crystal orientation which occurs in wetting layers at the freezing transition from the fluid 
phase at the walls (if complete wetting at the transition occurs). 

Thus, we have chosen values of D such that the thickness is compatible with an integer number 
of stacked (111) lattice planes without creating a noticeable elastic distortion of the crystal. Only 

11 



the thermodynamic integration method based on Eq. [T2]is used, and system hnear dimensions x 
Ly^D are taken, with L^xLy= [(8.8723 x 7.6835), (8.8331 x 7.6496), (8.7807 x 7.6043), (8.7290 x 
7.5595), (8.6779 x 7.5152), (8.6292 x 7.4730)] and several choices of D, corresponding to n = 
6, 12, 24, 48 stacked lattice planes. The result for does depend on D but is compatible with a 
linear variation in 1/D. So we find 7^c from an extrapolation versus 1/D ^ 0. We have checked 
the reliability of this approach for the case of the hard wall potential, Eq. [H where previous work 



with different methods have given [25| 7ic^(?7f) = 1.457 ± 0.018. 



Fig. [7] shows our results for 7^^^,^ (77) for the WCA potential as a function of packing fraction 
and several choices of e. The corresponding data for '-/^f for the fluid near the transition are also 
included. We find that the choice e = 1 yields functions 7^/(77), ■jwciv) which are very close to 



the corresponding data for the hard wall case. For the latter, Fortini and Dijkstra 25| have found 
that 'jwfivt) = 1-990 ± 0.007, and hence the difference ■jwfivt) ~ lli!'} = 0.53 ± 0.02. Laird and 
find 7^/(r/t) = 1.975 ± 0.002 and the difference 7^/(r/t) - 7^^,^ = 0.563 ± 0.004. 



Davidchack 



28 



These difference values are very close to the fiuid-crysta 
Laird and Davidchack using the cleaving method give [28 



interface tension. Recent simulations by 
7I11 = 0.557±0.007, 7^1° = 0.571±0.006 



and 7^^= 0.592 ± 0.007. The capillary wave fluctuation method, applied by the same authors 
gives 



29| 7I11 = 0.546 ± 0.016, 7"° = 0.557 ± 0.017 and 7^°° = 0.574 ± 0.017 and the most 
recent simulation resuhs from 2010 |4J] yield 7^^^ = 0.5416, 7^^° = 0.5590 and 7^°° = 0.5820 



with uncertainties in the last two digits. These results imply that complete wetting of the hard 
wall by the crystal in [111] orientation might occur {'jwfivt) — ll^ciVt) > 1^^^), but a finite, small 
contact angle cannot be excluded from the errorbars. We mention in passing that the difficulties 
in extracting interface tension with reliable errorbars might be substantial: as an example we 
mention the values for the interfacial stiffness 7^*^° obtained in different simulations using the 



capillary wave fluctuation method. Laird and Davidchack [29[ obtain 7^°° = 0.44 ± 0.03 (using 
thin slabs) whereas Zykova-Timan et al. obtain 7^°*^ = 0.49 ± 0.02 (using thick slabs). These 
stiffnesses include the anisotropy of the interfacial tension in an amplified manner but apparently 
depend on the simulation geometry. 

From Fig. [7] we conclude that changing the wall potential from the hard wall case {Eq. [1]} to 
the WCA case {Eq. [2l[ has little effect on the wetting properties of the wall, since the difference 
Iwfijit) ~ ll^ciVt) is independent of e, at least within the statistical errors of our calculation, and 
moreover is almost identical to the HW results. Thus, although the variation of e from e = 0.25 to 
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e = 4.0 enhances 'jwfivt) by about 0.5, the increase of 'jwfivt) is almost identical to the increase 
of 'Jwcivt), and hence one cannot reach a wetting transition (and then vary the contact angle) by 
varying e. 

VI. CONCLUSIONS 

In this work, the effects of confining walls on a hard sphere fluid were studied over a wide range 
of packing fractions, including also the regime of the transition to the solid crystalline phase. The 
effect of the wall was described by using a WCA potential {Eq. [2]} acting on the fluid particles, 
but for comparison also a hard wall potential {Eq. [T][ was chosen. The main interest of this 
paper was a comparative study of various methods to obtain the surface excess free energy and 
the surface excess density, applying both Monte Carlo (MC) methods and DFT calculations. We 
found very good agreement between all approaches in the fluid phase for not too large packing 
fractions {t] < 0.4), irrespective of the choice of the wall-fluid potential that was used. For 
7] > 0.35, systematic discrepancies between the MC and DFT results for the surface excess density 
were found, which presumably should be attributed to the fact that for high densities in the 
fluid nontrivial correlations between the fluid particles beyond the nearest neighbor shell develop, 
which are no longer described by DFT with very high accuracy. However, DFT describes very 
accurately the density distribution very close to the walls, and since this controls the wall-fluid 
surface tension, the latter is very accurately predicted by DFT (Fig. [5b). 

The application of MC methods for 77 > 0.4 also becomes increasingly difficult - the pronounced 
layering that occurs makes the procedures that we used sensitive to finite size effects both with 
respect to D (the regions disturbed by both walls start to interact) and with respect to L (when a 
precursor of a crystalline wetting layer occurs at a wall, the crystalline planes exhibit an in-plane 
triangular lattice structure, which exhibits a mismatch with an L x L cross-section due to the 
periodic boundary conditions). This problem occurs a fortiori in the solid phase (where also D 
needs to be chosen such that elastic distortion of the crystal in z-direction is avoided). Thus, 
our study is clearly a feasibility study only, and more work will be required to ascertain the true 
behavior occurring in the thermodynamic limit. We recall that for the case of hard walls, that we 
have included for comparison, many studies by different methods were indispensable to reach the 
current level of understanding. 
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One motivation of the present work was also to possibly control the difference 7iu/(e) — Iwd^) 
at the bulk fluid-solid transition by varying e, in order to allow a convenient study of a wetting 
transition at crystallization. However, unfortunately the variation of this difference with e is rather 
weak, and the system stays in the region of complete wetting (zero contact angle) or in the regime 
of small nonzero contact angles, so one cannot reach states deep in the incomplete wetting regime 
in this way. Nevertheless, our calculations could be useful to understand experiments were one 
uses walls coated with polymer brushes containing hard-sphere like colloidal dispersions. 

Acknowledgements: We acknowledge support by the Deutsche Forschungsgemeinschaft (DFG) 
under grants No Bi 314/19-2, SCHI 853/2-2, and SFB TR6 and the JSC for a grant of computer 
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Appendix A: Determination of surface tensions using density functional theory 

The equilibrium solvent density profile p(r) = Pcq(r) can be determined directly from the basic 
equations of density functional theory. The grand potential functional is given by 

m = ^'1p] + -^^^[P] - / drifi - V^'^'ir)) , (Al) 

where J-""^ and J-"*^^ denote the ideal and excess free energy functionals of the solvent. The chemical 
potential in the hard sphere fluid is denoted by /i and and the wall (hard or soft) defines the external 
potential V^°^* (given by Eqs. [H [2]) and which depends only on the Cartesian coordinate z. The 
exact form of the ideal part of the free energy is given by 

/3jrid[p] = J dr(3f'^ir) = J Sr p{r) {\n[p{r)k^] - l) . (A2) 

Here, A is the de-Broglie wavelength and /3 = 1/ {k-BT) is the inverse temperature. The equilibrium 
density profile Pcq(r) for the solvent at chemical potential /i = ln(pj, A^) -|- (corresponding 
to the bulk density pb) is found by minimizing the grand potential in Eq. (lAip : 

In + /3\/-t(z) = + . (A3) 

Ps dp{z) 

For an explicit solution, it is necessary to specify the excess part of the free energy. We employ 
fundamental measure functionals which represent the most precise functionals for the hard sphere 
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fluid. Specifically we employ: 



= J dv r{{n[piv)]}) , (A4) 
/3/-({n[p(r)]}) = -no ln(l - ^3) + ^1(^3) ''^'''"'"^"''^ + 

1-^3 

nl - 3n2 n2 ■ n2 + q;t I (n2 ■ nx ■ n2 - Tm|) 



</22(%) 

<^2 



247r(l - nsY 
2n3-ni + 2(1-^3)111(1-723) 

3^3 

^ 2^3-3^1 + 2^^ + 2(1-723)2 111(1-^3) 

3nl 



Here, /'^^ is a free energy density which is a function of a set of weighted densities {n(r)} = 
{no, ni,n2, n^, rii, 112, t^t} with four scalar, two vector and one tensorial weighted densities. These 
are related to the density profile p(r) by nair) = J dr'p{r') w°'{r — r'). The weight functions, 
{w(r)} = {w^ , , w"^ , , wt} , depend on the hard sphere radius R = a/2 as follows: 

w^ = e{R-\r\), w"" = SiR - \r\) , = ^ , w° - 



= yfiR - |r|) , = £^ ^ K).- = - I) - Ir|) • (A5) 

Setting = in Eq. (lA4p corresponds to neglecting the tensorial weighted density. This is the 
White Bear II (WBII) functional derived in Ref. 46j. This functional is most consistent with 
restrictions imposed by morphological thermodynamics {47 1, see below for a discussion what this 
means for the hard wall surface tension. Setting = 1 corresponds to the tensor modification 
(originally introduced in Ref. 3]) of WBII (WBII-T) which facilitates the hard sphere crystal 
description. Coexistence densities, bulk crystal free energies, density anisotropies in the unit cell 
and vacancy concentrations are described very well using WBII-T [49|. 

From the equilibrium density profiles pcq{z), the surface tension can be determined as the excess 
over bulk grand potential: 



POO 

7K] = / dz [fMz)] + r[p.^{z)] -{p- y-*(z))peq(^) - u:,] 

J zn 



(A6) 



where zo denotes the location of the wall and the grand potential density in the bulk is given by 
the negative pressure, Uh = —p- Both the WBII and the WBII-T functional are consistent with 
the Carnahan-Starling equation for p. 
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In the case of a hard wall, the surface tension can be determined from a scaled particle argument 



46 



SOU as follows: 



7sp 



dm 



ln(l-r/fe) , r/fe(2 + 3r/b - 2r/2) 



{n}={n,} ^ vr(l-r/,)2 



+ 2!ll^_::iIL^2>lL . (A7) 



This surface tension is taken with respect to the wall position zq being at the physical wall and 
not at the surface of exclusion z'q = Zq + a /2 where the wall potential jumps from infinity to 
zero. Here, the derivative of /^^ has to be evaluated with the bulk values for the set of weighted 
densities: 723,6 = rjb, n2,h = Q/f^Vb, ni^b = 3/(7rcr2) ^^^^ ^ 6/(7rcr^) r]b, ni^b = n2,b = ?^t,6 = with 
rjb = cr^Ti/Gpb denoting the bulk packing fraction. For a consistent functional, both expressions 
for the surface tension IA6I and IA7I should agree. The WBII functional is very consistent in this 
respect, as illustrated in Tab. [J and the WBII-T functional is only slightly less consistent. For 
packing fractions larger than 0.45 (close to freezing) the inconsistency becomes noticeable, this is 
also where we observe the largest deviations from the simulation results. The analytical 7sp is 
still closest to the simulation results. 

For soft walls, no analytical result can be derived. One would extrapolate from the hard wall 
results that 7[pcq] from the WBII functional will give slightly better results than 7[pcq] from the 
WBII-T functional. This is indeed what we have observed in comparison to the simulations. 
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TABLE I. Comparison of surface tension 7sp vs. 7[pcq] of hard spheres against a hard wall for various 
bulk packing fractions up to freezing. 
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FIG. 1. Density profile p{z) vs. z, for a box of linear dimensions L = 12.41786, D = 25.61184, total 
particle number N = 2866, and five choices of e, as indicated. The upper left inset shows the first peak of 
p{z) close to the left wall, resolved on a much finer abscissa scale; the upper right inset shows the density 
in the central part of the box, resolved on a much larger ordinate scale, to show that for the different 
values of e essentially the same bulk density pi, in the center of the film is obtained. 
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L = 12.41786, L ,= 12.41786, L =25.61184 
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N = 3545 Tl = 0.470 E = 0.000 
N = 3545 T] = 0.470 E = 1.000 
N = 3696 Tl = 0.490 E = 0.000 
N = 3696 Tl = 0.490 E = 1.000 
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FIG. 2. Same as Fig. 1, but for = 3545 and A'^ = 3696, respectively. Note that in both cases two 
choices of e are shown, namely e = (hard wall system) and e = 1.0, but on the scales of the plot 
these data coincide. Insert shows p{z) in the center of the film on magnified scales, to show that at the 
densities chosen in this figure the distance D chosen here is not large enough to render the two walls 
strictly noninteracting (the systematic density oscillations do not completely die out in the center of the 
film). 
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FIG. 3. Plot of the surface excess packing fraction {—r]s) versus the packing fraction ijh = phir/G, in the 
bulk for several choices of the strength e of the WCA potential due to the wall (Eq. 2). Some data are 
obtained from the same geometry as in Figs. [H [2l performing runs in the NVT ensemble (total density 
p and corresponding packing fraction rj being held constant). Data using other choices of L and D are 
included, to check for finite size effects. Triangles pointing to the right correspond to a geometry of L = 5 
and D = 40. Triangles pointing down correspond to L = 13 and D = 50. All the other symbols correspond 
to a geometry of L = 12.418 and D = 25.612. For comparison, also data for hard wall boundaries (Eq. 1) 
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30(1 for the 



are included, both from the present work and from the data of Laird and Davidchack 
excess volume vn which can be related to the surface excess packing fraction as —rjs = pbVNT^ Symbols 
are Monte Carlo data, and lines show the the corresponding DFT results. Here, full curves correspond 
to the White Bear II functional and broken curves correspond to the White Bear II (Tensor) functional. 
The functionals differ in their applicability to describe the fluid-crystal transition (see appendix). 
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FIG. 4. Wall potential Vwca{z), density p{z) and product Vy/cA{z)p{z) plotted vs. z, in the regime 
0.36 <z< 0.63, for both e = 0.1 and e = 1.0, for the case % = 0.42476 (a) and % = 0.46443 (b). Density 
functional results (using the White Bear II functional) 

are included, as full curves. 
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FIG. 5. Wall-fluid surface tension of the hard sphere fluid confined by hard walls plotted vs. packing 
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30| and present results, due to the use of 



fraction r/^ in the bulk. Symbols show literature data 
Eq. E] and the thermodynamic integration method based on Eq. W2[ respectively; lines show the result 
of our DFT calculation (full lines - White Bear II functional, broken lines ~ White Bear II (Tensor) 
functional). Note that a factor a^/ksT is put equal to unity in this figure and the following figures 
throughout. 
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FIG. 6. Wall-fluid surface tension 7^,^ plotted vs. packing fraction r/t, for the WCA wall potential 
{Eq. [S]}-, varying its strength e from e = 0.25 to e = 4, as indicated. Symbols are MC data, lines show 
the result of our DFT calculation (full lines - White Bear II functional, broken lines - White Bear II 
(Tensor) functional). 
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FIG. 7. Wall-fluid surface tension and wall-crystal surface tension 7^c plotted vs. packing fraction, 
for the WCA potential {Eq. [2]}, for conditions near the fluid-solid transition. As in Fig. [6l the strength 
e of the WCA potential is varied: e = 0.25, 0.5, 1.0, 2.0, and 4.0 (from bottom to top). The symbols are 
MC data obtained from the thermodynamic integration method based on Eq. error bars are from 
the linear fit as function of . At = 0.4896 the smallest system size was excluded from the linear 
extrapolation since a crystalline layer was visible at the wall when looking at snapshots. 
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